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Abstract. The aim of this work is to provide fast and accurate approx- 
imation schemes for the Monte-Carlo pricing of derivatives in the Levy 
LIBOR model of Eberlein and Ozkan (2005) . Standard methods can be 
applied to solve the stochastic differential equations of the successive 
LIBOR rates but the methods are generally slow. We propose an alter- 
native approximation scheme based on Picard iterations. Our approach 
is similar in accuracy to the full numerical solution, but with the feature 
that each rate is, unlike the standard method, evolved independently of 
the other rates in the term structure. This enables simultaneous calcula- 
tion of derivative prices of different maturities using parallel computing. 
We include numerical illustrations of the accuracy and speed of our 
method pricing caplets. 

1. Introduction 

The LIBOR market model has become a standard model for the pric- 
ing of interest rate derivatives in recent years. The main advantage of the 
LIBOR model in comparison to other approaches, is that the evolution of 
discretely compounded, market-observable forward rates is modeled directly 
and not deduced from the evolution of unobservable factors. Moreover, the 
log-normal LIBOR model is consistent with the market practice of pricing 
caps according to Black's formula (cf. Black 1976). However, despite its ap- 
parent popularity, the LIBOR market model has certain well-known pitfalls. 

On the one hand, the log-normal LIBOR model is driven by a Brownian 
motion, hence it cannot be calibrated adequately to the observed market 
data. An interest rate model is typically calibrated to the implied volatility 
surface from the cap market and the correlation structure of at-the-money 
swaptions. Several extensions of the LIBOR model have been proposed in the 
literature using jump-diffusions. Levy processes or general semimartingales 
as the driving motion (cf. Glasserman and Kou 2003, Eberlein and Ozkan 
2005, Jamshidian 1999), or incorporating stochastic volatility effects (cf. e.g. 
Andersen and Brotherton-Ratcliffe 2005). 

On the other hand, the dynamics of LIBOR rates are not tractable un- 
der every forward measure due to the random terms that enter the dy- 
namics of LIBOR rates during the construction of the model. In particular, 
when the driving process has continuous paths the dynamics of LIBOR rates 
are tractable under their corresponding forward measure, but they are not 
tractable under any other forward measure. When the driving process is 
a general semimartingale, then the dynamics of LIBOR rates are not even 
tractable under their very own forward measure. Consequently: if the driving 
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process is a continuous semimartingale caplets can be priced in closed form, 
but not swaptions or other multi-LIBOR derivatives. However, if the driving 
process is a general semimartingale, then even caplets cannot be priced in 
closed form. The standard remedy to this problem is the so-called "frozen 
drift" approximation, where one replaces the random terms in the dynamics 
of LIBOR rates by their deterministic initial values; it was first proposed by 
Brace et al. (1997) for the pricing of swaptions and has been used by several 
authors ever since. Brace et al. (2001) among others argue that freezing the 
drift is justified, since the deviation from the original equation is small in 
several measures. 

Although the frozen drift approximation is the simplest and most popular 
solution, it is well-known that it does not yield acceptable results, especially 
for exotic derivatives and longer horizons. Therefore, several other approxi- 
mations have been developed in the literature. We refer the reader to Joshi 
and Stacey (2008) for a detailed overview of that literature, and for some 
new approximation schemes and numerical experiments. 

Although most of this literature focuses on the lognormal LIBOR market 
model, Glasserman and Merener (2003b, 2003a) have developed approxima- 
tion schemes for the pricing of caps and swaptions in jump-diffusion LIBOR 
market models. 

In this article we develop a general method for the approximation of the 
random terms that enter into the drift of LIBOR models. In particular, 
by applying Picard iterations we develop a generic approximation scheme. 
The method we develop yields more accurate results than the frozen drift 
approximation, while having the added feature that the individual rates can 
be evolved independently in a Monte Carlo simulation. This enables the 
use of parallel computing in the maturity dimension. Moreover, our method 
is universal and can be applied to any LIBOR model driven by a general 
semimartingale. We illustrate the accuracy and speed of our method in a 
case where LIBOR rates are driven by a normal inverse Gaussian process. 



2. The Levy LIBOR model 

The Levy LIBOR model was developed by Eberlein and Ozkan (2005), 
following the seminal articles of Sandmann et al. (1995), Miltersen et al. 
(1997) and Brace et al. (1997) on LIBOR market models driven by Brow- 
nian motion; see also Glasserman and Kou (2003) and Jamshidian (1999) 
for LIBOR models driven by jump processes and general semimartingales 
respectively. The Levy LIBOR model is a market model where the forward 
LIBOR rate is modeled directly, and is driven by a time-inhomogeneous 
Levy process. 

Let = Tq < Ti < • • ■ < T/v < T^+i = denote a discrete tenor 
structure where 6i = Tj+i — Ti, i £ {0,1, ... , N}. Consider a complete sto- 
chastic basis (ri, J^, FjPt^^) and a time-inhomogeneous Levy process H = 
{Ht)o<t<T, satisfying standard assumptions such as the existence of expo- 
nential moments and absolutely continuous characteristics. The law of H is 
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described by the Levy-Khintchine formula: 

t 

JEpt, [e'"^'] = exp y Ks{iu)d^ . (2.1) 



Here is the cumulant generating function associated to the infinitely di- 
visible distribution with Levy triplet (0, c, F'^*), i.e. for u G M and s G [0, T*] 

K,{iu) = -^u^ + J (e™^ - 1 - iux)F^'{dx). (2.2) 
The canonical decomposition of H is: 

H = j ^dVFj'+ j j x(/i^ -z^^*)(ds,dx), (2.3) 

OR 

where W'^* is a P^, -standard Brownian motion, /i^ is the random measure 
associated with the jumps of H and v"^* is the Pt, -compensator of //^. We 
further assume that the following conditions are in force. 

(LRl): For any maturity Tj there exists a bounded, continuous, deter- 
ministic function A(-,Tj) : [0, Tj] — ?• M, which represents the volatil- 
ity of the forward LIBOR rate process L(-,Tj). Moreover, we as- 
sume that (i) for all s G [0,T*], there exist M, e > such that 
/{|x|>i}e'''^^t(da;)dt < oo, for u G [-(l + e)M, (l + e)M], and (ii) 
for all s < Ti 



N 



.Ti) < M. 



i=l 



(LR2): The initial term structure B{0,Ti), 1 < i < + 1, is strictly 
positive and strictly decreasing. Consequently, the initial term struc- 
ture of forward LIBOR rates is given, for 1 < i < A^, by 

The construction of the model starts by postulating that the dynamics 
of the forward LIBOR rate with the longest maturity L{-,Tn) is driven by 
the time-inhomogeneous Levy process H and evolve as a martingale under 
the terminal forward measure Pt, • Then, the dynamics of the LIBOR rates 
for the preceding maturities are constructed by backward induction; they 
are driven by the same process H and evolve as martingales under their 
associated forward measures. For the full mathematical construction we refer 
to Eberlein and Ozkan (2005). 

We will now proceed to introduce the stochastic differential equation that 
the dynamics of log-LIBOR rates satisfy under the terminal measure Pt* . 
This will be the starting point for the approximation method that will be 
developed in the next section. 
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In the Levy LIBOR model the dynamics of the LIBOR rate L(-, Tj) under 
the terminal forward measure JPt, are given by 

L(t,T,) = L(0,ri)exp b{s,T^)ds + J A(s,r,)di?,j , (2.5) 

where H = {Ht)o<t<T, is the Pr^-time-inhomogeneous Levy process. The 
drift term b{-,Ti) is determined by no-arbitrage conditions and has the form 

b{s,T,) = -lx\s,T,)cs-CsX{s,T,) ^ -A^pI^X(s,Ti) 



gA(s,T,)x 



TV \ 

l) n (3{s,x,Ti)-Xis,Ti)x\ F^'idx), (2.6) 

l=i+l J 



where 

Note that the drift term in (2.5) is random, therefore we are dealing with a 
general semimartingale, and not with a Levy process. Of course, L(-,Tj) is 
not a P^^, -martingale, unless i = N (where we use the conventions "^1=1 = 

and Uli = !)• 

Let us denote by Z the log-LIBOR rates, that is 

Z{t,T,) ■.= \og L{t,Ti) 

t t 
= Z{0,Ti) + J bis,Ti)ds + J Xis,Ti)dHs, (2.8) 

where Z{0,Ti) =logL(0,Ti) for alH G {!,..., A^}. 

Remark 2.1. Note that the martingale part of Z(-,Tj), i.e. the stochastic 
integral X{s,Ti)dHs, is a time-inhomogeneous Levy process. However, the 
random drift term destroys the Levy property of Z(-,Tj), as the increments 
are no longer independent. 



3. PiCARD APPROXIMATION FOR LIBOR MODELS 

The log-LIBOR can be alternatively described as a solution to the follow- 
ing linear SDE 

dZ(t, Ti) = bit, Ti)dt + X{t, T,)dHt, (3.1) 

with initial condition Z(0, Tj) = logL(0,Tj). Let us look further into the 
above SDE for the log-LIBOR rates. We introduce the term Z(-) in the drift 
term b{-,Ti; Z(-)) to make explicit that the log-LIBOR rates depend on all 
subsequent rates on the tenor. 

The idea behind the Picard approximation scheme is to approximate 
the drift term in the dynamics of the LIBOR rates; this approximation 
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is achieved by the Picard iterations for (3.1). The first Picard iteration for 
(3.1) is simply the initial value, i.e. 

z(o)(t,ri) = z(o,ri), (3.2) 

while the second Picard iteration is 

t t 

Z^^\t,T,) = Z{0,T,) + J b{s,Tf,Z^'^\s))ds + J X{s,Ti)dH, 





t t 



= Z(0, T,) + J b{s, Ti- Z(0))ds + j \{s, Ti)dHs. (3.3) 



Since the drift term b{-,Ti; Z(0)) is deterministic, as the random terms have 
been replaced with their initial values, we can easily deduce that the second 
Picard iterate Z^^\-,Ti) is a Levy process. 

Comparing (3.3) with (2.8) it becomes evident that we are approximat- 
ing the semimartingale Z(-,Tj) with the time-inhomogeneous Levy process 

3.1. Application to LIBOR models. In this section, we will apply the Pi- 
card approximation of the log-LIBOR rates Z(-, Tj) by Z^^\-,Ti) in order to 
derive a strong, i.e. pathwise, approximation for the dynamics of log-LIBOR 
rates. That is, we replace the random terms in the drift 6(-,Tj; Z{-)) by the 
Levy process Z^^\-,Ti) instead of the semimartingale Z(-,Tj). Therefore, 
the dynamics of the approximate log-LIBOR rates are given by 

t t 

Z{t,T,) = ZiO,Ti) + J b{s,Ti;Z^^\s))ds + J Xis,Ti)dHs, (3.4) 



where the drift term is provided by 

b{s,Tf,Z^^Hs)) = --X\s,Ti)c, - CsXis,T,) ^ \ ^, X{s,Ti) 

N 

X{s,Ti)x 




N \ 

l) J] d{s,x,Ti)-X{s,T,)x\F^'{da 

l=i+l J 



(3.5) 

with 

^TTi^S§^(«^"'^'"-)-- 

The main advantage of the Picard approximation is that the resulting 
SDE for Z{-,Ti) can be simulated more easily than the equation for Z(-, Tj). 
Indeed, looking at (3.1) and (2.6) again, we can observe that each LIBOR 
rate L(-, Tj) depends on all subsequent rates L(-, Ti), i + 1 < I < N . Hence, 
in order to simulate L{-,Ti), we should start by simulating the furthest rate 
in the tenor and proceed iteratively from the end. On the contrary, the 
dynamics of Z(-,Tj) depend only on the Levy processes Z^^\-,Ti), i + 1 < 
I < N, which are independent of each other. Hence, we can use parallel 
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computing to simulate all approximate LIBOR rates simultaneously. This 
significantly increases the speed of the Monte Carlo simulations which will 
be demonstrated in the numerical example. 

3.2. Drift expansion. Let us look now at the drift term in (2.6) more 
carefully. Observe that there is a product of the form nfc=i(-'-+'^fc) appearing; 
the expansion of this product has the following form 

N N 

JJ (1 + Ofc) = 1 + ^ Ofc + ^ aiUj + ^ aittjak 

k=l k=l l<i<j<N l<i<j<k<N 

N 

+ ^ aiaja^ai H h 

i^j^k^l k=l 
' V ' 

where the number of terms on the right hand side is 2^. Therefore, we need 
to perform 2^ computations in order to calculate the drift of the LIBOR 
rates. Since N is the length of the tenor, it becomes apparent that this 
calculation is feasible for a short tenor, but not for long tenors; e.g. for 

= 40 this amounts to more than 1 trillion computations. 

In order to deal with this computational problem, we will approximate 
the LHS of (3.7) with the first or second order terms. Let us introduce the 
following shorthand notation for convenience: 

Xr.= X{s,Ti) and Li := L{s,Ti). (3.8) 

We denote by A the part of the drift term that is stemming from the jumps, 
i.e. 

N 

A 



/(p._l)n_(l + _^(e.._,))_A..]Fj-(d.). 

(3.9) 



The first order approximation of the product term is 



A'=l{ (e^- - l) (^1 + E^-^y- - l) ) - A.. ) (dx) 



e^«^ - 1 - Xixj fJ* (dx) 



1=1+1 

N 



i=i+i '■ 
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and the order of the error is 

A = A' + 0(||L|p). (3.11) 
Similarly the second order approximation is provided by 

i=i+i ' '■ 

i+l<k<l<N 

K{Xi + Xl + Afc) - K{Xi + A;) - K{Xi + Xk) 

- K{Xk + Xi) + K{Xi) + k{Xi) + K(Afc)) , (3.12) 
and the order of the error is 

A = A" + 0(||Lf ). (3.13) 

3.3. Caplets. The price of a caplet with strike K maturing at time Tj, 
using the relationship between the terminal and the forward measures can 
be expressed as 

N 

CoiK,Ti) = 5B{0,T,)lEp^^[ J] {l + 6L{Ti,Ti)){L{Ti,Ti) - K)^ . 

l=i+l 

(3.14) 

This equation will provide the actual prices of caplets corresponding to sim- 
ulating the full SDE for the LIBOR rates. In order to calculate the Picard 
approximation prices for a caplet we have to replace L(-,T.) in (3.14) with 
L{-,T.). Similarly, for the frozen drift approximation prices we must use 
L°(-,T.) instead of L{-,T.). 

4. Numerical illustration 

The aim of this section is to demonstrate the accuracy and efficiency of 
the Picard approximation scheme for the valuation of options in the Levy 
LIBOR model compared to the "frozen drift" approximation. In addition, 
we investigate the accuracy of the drift expansions in section 3.2. We will 
consider the pricing of caplets, although many other interest rate derivatives 
can be considered in this framework. 

We revisit the numerical example in Kluge (2005, pp. 76-83). That is, we 
consider a tenor structure Tq = 0, Ti = ^,T2 = I . . . , Tiq = 5 = T^,, constant 
volatilities 

A(-,ri) = 0.20 A(-,r2) = 0.19 A(-,r3) = 0.18 
A(-, T4) = 0.17 A(-, T5) = 0.16 A(-, Te) = 0.15 
X{-,Tj) = 0.14 A(-, Ts) = 0.13 A(-, Tg) = 0.12 

and the discount factors (zero coupon bond prices) as quoted on February 
19, 2002; cf. Table 4.1. The tenor length is constant and denoted by (5 = |. 

The driving Levy process H is a normal inverse Gaussian (NIG) process 
with parameters a = 6 = 1.5 and fi = P = 0. We denote by fi^ the random 
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T 


0.5Y 


lY 


1.5Y 


2Y 


2.5Y 


BiO,T) 


0.9833630 


0.9647388 


0.9435826 


0.9228903 


0.9006922 


T 


3Y 


3.5Y 


4Y 


4.5Y 


5Y 


i?(0,T) 


0.8790279 


0.8568412 


0.8352144 


0.8133497 


0.7920573 



Table 4.1. Euro zero coupon bond prices on February 19, 2002. 



measure of jumps of H and by iy{dt,dx) = F{dx)dt the Pt» -compensator 
of /i^, where F is the Levy measure of the NIG process. The necessary 
conditions are satisfied because M = a, hence | A(-, T!j)| = 1.44 < a 

and A(-, Ti) < f , for alH G {1, . . . , 9}. 

The NIG Levy process is a pure-jump Levy process with canonical de- 
composition H = /q /jj x{fi^ — u){ds, dx). The cumulant generating function 
of the NIG distribution, for all u G C with \?R.u\ < a, is 

k{u) = 6a — b\l o? — V?. (4.1) 

In figure 4.1 we plot the difference in basis points (bp) between caplet 
implied volatility calculated from the full numerical solution and implied 
volatilities from the frozen drift and the Picard approximation respectively. 
In order to isolate the error from the two approximations we use the same 
discretization grid (5 steps per tenor length) and the same pseudo random 
numbers (10000 paths) in each method. The pseudo random numbers are 
generated from the NIG distribution using the standard methodology de- 
scribed in Glasserman 2003. The drifts are first calculated without approx- 
imation using expression (3.5). It is clear that the Picard approximation 
outperforms the frozen drift with an error which is at maximum 0.023 bp. 
Since implied volatility is quoted in units of a basis point then Ibp is a 
natural maximum tolerance level of error in an approximation. It is clear 
that the error in the frozen drift approximation is significantly bigger than 
one basis point throughout most strikes and maturities, with a maximum 
around 17bp. Both graphs also show that the error in general increases in 
absolute terms the smaller the strike. 

As we established in (3.7) the number of terms needed to calculate the 
drift grows with a rate 2^. In market applications is often as high as 
60 reflecting a 30 year term structure with a 6 month tenor increment. At 
this level even the calculation of one drift term becomes infeasible and this 
necessitates the use of the approximations introduced in (3.10) and (3.12). 
If we investigate the errors introduced by comparing them with the full 
numerical solution we get an average error of 0.41 bp with max of 9.5 whereas 
the second order approximation performs much better with an average error 
of 0.013 and maximum error around 0.38 bp. 

In terms of computational time a large gain is realized when using the 
approximations in (3.10) and (3.12). In the example above the CPU time 
for the full numerical solution is 141 seconds but after applying the first or- 
der or second order drift approximation it drops to 0.9 and 1.2 respectively. 
Adding the Picard approximation to these three cases does not contribute 
to the computational speed unless parallelization is employed. On the left 
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Figure 4.1. Difference in implied caplet volatility (in basis 
points) between the full SDE and the frozen drift prices (left), 
and the full SDE and the Picard approximation (right). 




Figure 4.2. CPU time as function of the number of paths 
(left) and CPU time as a function of N, the number of rates 
(right). 



in figure 4.2, CPU time as function of the number of paths for the Picard 
approximation and the full numerical solution is plotted. Both use the sec- 
ond order drift approximation scheme in (3.12). The computations are done 
in Matlab running on an Intel 17 processor with the capability of running 
8 processes simultaneously. Here we see the typical linear behavior as the 
number of paths are increased but it can be seen that the Picard approxi- 
mation has a significantly lower slope. Furthermore, on the right when we 
plot CPU time as a function of rates one can see CPU time exponentially 
increasing, revealing that large gains in computational time are realizable 
when using the Picard approximation scheme and the drift expansion. 

5. Conclusion 

This paper derives a new approximation method for Monte Carlo deriv- 
ative pricing in LIBOR models. It is generic and can be used for any semi- 
martingale driven model. It decouples the interdependence of the rates when 
moving them forward in time in a simulation, meaning that the computa- 
tions can be parallelized in the maturity dimension. We have demonstrated 
both the accuracy and speed of the method in a numerical example. 
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